// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtReflectFilter.C
// ************************************************************************* //

#include <avtReflectFilter.h>

#include <float.h>
#include <VisItInit.h>
#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkIdList.h>
#include <vtkPointData.h>
#include <vtkPointSet.h>
#include <vtkRectilinearGrid.h>
#include <vtkUnsignedCharArray.h>

#include <avtDataTree.h>
#include <avtDatasetExaminer.h>
#include <avtExtents.h>

#include <DebugStream.h>
#include <ImproperUseException.h>

#include <vtkStructuredGrid.h>
#include <vtkRectilinearGridToPointSet.h>

using   std::string;

static void ReflectVectorData(vtkDataSet *in_ds, int dim, bool);

// ****************************************************************************
//  Function: Equal
//
//  Purpose:
//      Compares a float and a double and determines if they are equal.
//
//  Note:   Tried to do this with integer arithmetic but ran into compiler bugs
//          with gcc2.95 (worked with print statements, not without, with -O2).
//
//  Programmer: Hank Childs
//  Creation:   August 11, 2004
//
// ****************************************************************************

static inline bool
Equal(float t1, double t2)
{
    float v1 = t1;
    float v2 = (float) t2;

    // If they are exactly even, then no need to do more.
    if (v1 == v2)
        return true;

    // Rule out the 0's.
    if (v1 == 0.)
        return (fabs(v2) < 1e-10 ? true : false);
    if (v2 == 0.)
        return (fabs(v1) < 1e-10 ? true : false);
 
    // Make sure they are signed correctly.
    if (v1 < 0. && v2 > 0.)
        return false;
    if (v1 > 0. && v2 < 0.)
        return false;

    v1 = fabs(v1);
    v2 = fabs(v2);

    float greater = (v1 > v2 ? v1 : v2);
    float diff = fabs(v1-v2);
   
    // Save ourselves a divide -- below is the same as diff/greater < 1e-6.
    if (diff < (1e-6)*greater)
        return true;

    return false;
}


// ****************************************************************************
//  Method: avtReflectFilter constructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Mar 7 10:35:24 PDT 2002
//
//  Modifications:
//
//    Hank Childs, Fri Sep  3 12:10:47 PDT 2010
//    Add Boolean for zero-ing velocities on the boundary.
//
//    Alister Maguire, Tue Apr 17 11:12:30 PDT 2018
//    Added ghostsCreated. 
//
// ****************************************************************************

avtReflectFilter::avtReflectFilter()
{
    xReflect = 0.;
    yReflect = 0.;
    zReflect = 0.;
    zeroOutVelocitiesOnBoundary = false;
    ghostsCreated = false;
}


// ****************************************************************************
//  Method: avtReflectFilter destructor
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Mar 7 10:35:24 PDT 2002
//
//  Modifications:
//
// ****************************************************************************

avtReflectFilter::~avtReflectFilter()
{
}


// ****************************************************************************
//  Method:  avtReflectFilter::Create
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Mar 7 10:35:24 PDT 2002
//
// ****************************************************************************

avtFilter *
avtReflectFilter::Create()
{
    return new avtReflectFilter();
}


// ****************************************************************************
//  Method:      avtReflectFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Mar 7 10:35:24 PDT 2002
//
// ****************************************************************************

void
avtReflectFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const ReflectAttributes*)a;
}


// ****************************************************************************
//  Method: avtReflectFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtReflectFilter with the given
//      parameters would result in an equivalent avtReflectFilter.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Mar 7 10:35:24 PDT 2002
//
// ****************************************************************************

bool
avtReflectFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(ReflectAttributes*)a);
}


// ****************************************************************************
//  Method: avtReflectFilter::PreExecute
//
//  Purpose:
//      Determines what lines/planes the data should be reflected across.  This
//      may be trivial (if they are all specified) or require some work (if
//      we are running in parallel).
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2002
//
//  Modifications:
//
//    Jeremy Meredith, Mon Mar 11 15:41:55 PST 2002
//    Made octant 0 be +x+y+z.
//
//    Hank Childs, Mon Feb 24 17:19:42 PST 2003
//    Better checking of when we need to calculate extents in parallel.
//
//    Hank Childs, Wed Jun 25 09:30:59 PDT 2003
//    Fix stupid cut-n-paste bug.
//
//    Jeremy Meredith, Thu Feb 15 11:55:03 EST 2007
//    Call inherited PreExecute before everything else.
//
// ****************************************************************************

void
avtReflectFilter::PreExecute(void)
{
    avtSIMODataTreeIterator::PreExecute();

    //
    // Calculating the extents may require parallel communication.  Only do
    // this if we are sure that we will need the extents.
    //
    double extents[6] = { -DBL_MAX, DBL_MAX, -DBL_MAX, DBL_MAX, 
                          -DBL_MAX, DBL_MAX };
    bool reflectX = false;
    bool reflectY = false;
    bool reflectZ = false;
    const int *shouldreflect = atts.GetReflections();
    for (int i = 0 ; i < 8 ; i++)
    {
        if (shouldreflect[i] != 0)
        {
            if (i & 1)
                reflectX = true;
            if (i & 2)
                reflectY = true;
            if (i & 4)
                reflectZ = true;
        }
    }
    bool needX = reflectX && atts.GetUseXBoundary();
    bool needY = reflectY && atts.GetUseYBoundary();
    bool needZ = reflectZ && atts.GetUseZBoundary();
    if (needX || needY || needZ)
    {
        GetSpatialExtents(extents);
    }

    if (atts.GetUseXBoundary())
    {
        //
        // If we are in an octant that is positive for X, then we want to
        // reflect along the minimum x-value.  Also, the octant number is set
        // up such that every other octant is positive for X.
        //
        int octant = atts.GetOctant();
        xReflect = (octant & 1 ? extents[1] : extents[0]);
    }
    else
    {
        xReflect = atts.GetSpecifiedX();
    }

    if (atts.GetUseYBoundary())
    {
        //
        // If we are in an octant that is positive for Y, then we want to
        // reflect along the minimum y-value.  Also, the octant number is set
        // up such that positive Y is paired off in twos.
        //
        int octant = atts.GetOctant();
        yReflect = (octant & 2 ? extents[3] : extents[2]);
    }
    else
    {
        yReflect = atts.GetSpecifiedY();
    }

    if (atts.GetUseZBoundary())
    {
        //
        // If we are in an octant that is positive for Y, then we want to
        // reflect along the minimum y-value.  Also, the octant number is set
        // up such that positive Z is the highest 4.
        //
        int octant = atts.GetOctant();
        zReflect = (octant & 4 ? extents[5] : extents[4]);
    }
    else
    {
        zReflect = atts.GetSpecifiedZ();
    }
}


// ****************************************************************************
//  Method: avtReflectFilter::PostExecute
//
//  Purpose:
//      Cleans up after the execution.  This manages extents.
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2002
//
//  Modifications:
//
//    Hank Childs, Fri Mar 15 19:49:31 PST 2002
//    Account for dataset examiner.
//
//    Hank Childs, Thu Feb 26 16:19:33 PST 2004
//    Fix typo where data extents are being reset.
//
//    Hank Childs, Wed Aug 11 15:09:07 PDT 2004
//    Tell the output that it has ghost zones.
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
//    Alister Maguire, Wed Apr 11 09:29:49 PDT 2018
//    Added check to make sure ghost zones have been created
//    before setting atts. 
//
// ****************************************************************************

void
avtReflectFilter::PostExecute(void)
{
    avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();
    outAtts.GetOriginalSpatialExtents()->Clear();
    outAtts.GetDesiredSpatialExtents()->Clear();
    outAtts.GetActualSpatialExtents()->Clear();

    if (GetInput()->GetInfo().GetAttributes().GetContainsGhostZones()
           != AVT_HAS_GHOSTS && ghostsCreated)
        outAtts.SetContainsGhostZones(AVT_CREATED_GHOSTS);

    double bounds[6];
    avtDataset_p ds = GetTypedOutput();
    avtDatasetExaminer::GetSpatialExtents(ds, bounds);
    outAtts.GetThisProcsOriginalSpatialExtents()->Set(bounds);
}


// ****************************************************************************
//  Method: avtReflectFilter::ModifyContract
//
//  Purpose:
//      Performs a restriction on the specification.  This is in place to
//      insure that we do static load balancing when extent calculation is
//      required.
//
//  Arguments:
//      spec    The current specification.
//
//  Returns:    A specification with the balancing potentially restricted.
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2002
//
//  Modifications:
//    Kathleen Bonnell, Wed Jun 19 12:28:10 PDT 2002
//    Don't turn off Zone numbers if they have been turned on elsewhere in
//    the pipeline.
//
//    Kathleen Bonnell, Wed Jun 19 13:42:37 PDT 2002
//    Completely removed the code turning off zone numbers.  Why set a flag
//    to false if it is already false?  False is the default setting.
//
//    Kathleen Bonnell, Wed Jun  2 09:21:46 PDT 2004 
//    Turn on node numbers when required. 
//
//    Hank Childs, Tue Feb 19 19:45:43 PST 2008
//    Rename "dynamic" to "streaming", since we really care about whether we
//    are streaming, not about whether we are doing dynamic load balancing.
//    And the two are no longer synonymous.
//
//    Hank Childs, Fri Sep  3 12:10:47 PDT 2010
//    Set Boolean for zero-ing velocities on the boundary.
//
// ****************************************************************************

avtContract_p
avtReflectFilter::ModifyContract(avtContract_p spec)
{
    avtContract_p ns = new avtContract(spec);

    //
    // If we need to calculate the extents to reflect across, then we will have
    // to have all the domains together at the same time, so we have to use
    // static load balancing.
    //
    if (atts.GetUseXBoundary() || atts.GetUseYBoundary() ||
        atts.GetUseZBoundary())
    {
        ns->NoStreaming();
    }

    if (spec->GetDataRequest()->MayRequireZones())
    {
        ns->GetDataRequest()->TurnZoneNumbersOn();
    }
    if (spec->GetDataRequest()->MayRequireNodes())
    {
        ns->GetDataRequest()->TurnNodeNumbersOn();
    }

    zeroOutVelocitiesOnBoundary = ns->GetDataRequest()->VelocityFieldMustBeContinuous();

    return ns;
}


// ****************************************************************************
//  Method: avtReflectFilter::UpdateDataObjectInfo
//
//  Purpose:
//     Specifies that this filter transforms the original dataset points. 
//
//  Programmer: Kathleen Bonnell 
//  Creation:   March 26, 2002 
//
//  Modifications:
//    Kathleen Bonnell, Mon Apr 14 09:57:39 PDT 2003
//    Set CanUseTransform to false.
//
//    Kathleen Bonnell, Wed Jun  2 09:21:46 PDT 2004 
//    Set CanUseInvTransform to false.
// 
//    Hank Childs, Fri Jan 13 09:58:47 PST 2006
//    Invalidate spatial meta-data.
//
//    Hank Childs, Tue Nov 30 11:56:13 PST 2010
//    Indicate that the zones have been invalidated.  If we do a rectilinear 
//    reflection, then we re-order the cells.  
//
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************
 
void
avtReflectFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().SetPointsWereTransformed(true);

    //
    // This filter invalidates any transform matrix in the pipeline.
    //
    GetOutput()->GetInfo().GetAttributes().SetCanUseInvTransform(false);
    GetOutput()->GetInfo().GetAttributes().SetCanUseTransform(false);

    // There are more zones and nodes than originally, so invalidate
    // the original information
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().InvalidateNodes();

    GetOutput()->GetInfo().GetValidity().InvalidateSpatialMetaData();

    GetOutput()->GetInfo().GetAttributes().AddFilterMetaData("Reflect");
}


// ****************************************************************************
//  Method: avtReflectFilter::ExecuteDataTree
//
//  Purpose:
//      Sends the specified input and output through the Reflect filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data tree.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Thu Mar 7 10:35:24 PDT 2002
//
//  Modifications:
//    Eric Brugger, Thu Jul 31 19:15:48 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Alister Maguire, Tue Apr 10 13:22:11 PDT 2018
//    Added the option to reflect across an arbitrary plane. 
//
// ****************************************************************************

avtDataTree_p 
avtReflectFilter::ExecuteDataTree(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set, the domain number, and the label.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();
    int domain = in_dr->GetDomain();
    std::string label = in_dr->GetLabel();

    if (in_dr == NULL)
    {
        return NULL;
    }

    avtDataTree_p rv = NULL;
    int reflectType  = atts.GetReflectType();

    switch(reflectType)
    {
        case ReflectAttributes::Plane:
        {
            //
            // For reflecting about an arbitrary plane, we
            // only need two output datasets: the original
            // and the reflected. 
            //
            vtkDataSet *reflections[2];
            reflections[0] = PlaneReflect(in_ds, false);
            reflections[1] = PlaneReflect(in_ds, true);
   
            //
            // Construct a data tree out of our reflections.
            //
            rv = new avtDataTree(2, reflections, domain, label);

            //
            // Clean up memory.
            //
            for (int i = 0 ; i < 2 ; i++)
            {
                if (reflections[i] != NULL)
                {
                    reflections[i]->Delete();
                }
            }
        }
        break;

        case ReflectAttributes::Axis:
        {
            //
            // Perform any of the 8 possible reflections.
            //
            int  i;
            const int *shouldreflect = atts.GetReflections();
            vtkDataSet *reflections[8];
            for (i = 0 ; i < 8 ; i++)
            {
                reflections[i] = NULL;
                if (shouldreflect[i] != 0)
                {
                    reflections[i] = AxisReflect(in_ds, i);
                }
            }

            //
            // Construct a data tree out of our reflections.
            //
            rv = new avtDataTree(8, reflections, domain, label);

            //
            // Clean up memory.
            //
            for (i = 0 ; i < 8 ; i++)
            {
                if (reflections[i] != NULL)
                {
                    reflections[i]->Delete();
                }
            }
        }
        break;
    }

    return rv;
}


// ****************************************************************************
//  Method: avtReflectFilter::PlaneReflect
//
//  Purpose:
//      If doReflection is true, create a new dataset from the input
//      and reflect this dataset across a given plane. Otherwise, 
//      create a new dataset that is identical to the old one but 
//      a vtkPointSet. 
//
//  Arguments:
//      in_ds          The input data. 
//      doReflection   Whether or not we should reflect the input data. 
//
//  Returns:    The dataset appropriately reflected.
//
//  Programmer: Alister Maguire
//  Creation:   Mon Apr  9 12:50:15 PDT 2018
//
//  Modifications:
//
// ****************************************************************************

vtkDataSet *
avtReflectFilter::PlaneReflect(vtkDataSet *in_ds, bool doReflection)
{
    vtkPointSet *out_ds = NULL;
    vtkRectilinearGridToPointSet *gridToPointFilter = NULL;
    
    //
    // If our input data is a rectilinear grid, we need to convert 
    // this to a point set. Otherwise, we can just cast the input
    // to a point set. 
    //
    if (in_ds->GetDataObjectType() == VTK_RECTILINEAR_GRID)
    {
        out_ds = (vtkPointSet *)vtkStructuredGrid::New();
        vtkRectilinearGridToPointSet *gridToPointFilter = 
            vtkRectilinearGridToPointSet::New();
        gridToPointFilter->SetInputData(in_ds);
        gridToPointFilter->Update();
        out_ds->ShallowCopy(gridToPointFilter->GetOutput());
    }
    else
    {
        out_ds = (vtkPointSet *)in_ds->NewInstance();
        out_ds->ShallowCopy(in_ds);
    }

    if (!doReflection)
        return out_ds;

    double *planePoint  = atts.GetPlanePoint();
    double *planeNormal = atts.GetPlaneNormal();
    double  beta        = 0.0;
    double  normalMag   = sqrt(planeNormal[0]*planeNormal[0]+
                               planeNormal[1]*planeNormal[1]+
                               planeNormal[2]*planeNormal[2]);

    //
    // Make sure the normal is normalized, and calculate
    // beta in the plane equation ax + by + cz = beta. 
    //
    for (int i = 0; i < 3; ++i)
    {
        planeNormal[i] = planeNormal[i] / normalMag;
        beta += planeNormal[i]*planePoint[i];
    }

    vtkPoints *inPts  = out_ds->GetPoints();
    vtkPoints *outPts = vtkPoints::New(inPts->GetDataType());
    int nPts          = inPts->GetNumberOfPoints();
    outPts->SetNumberOfPoints(nPts);

    for (int i = 0 ; i < nPts ; i++)
    {
        double reflection[3];
        double curPt[3];
        inPts->GetPoint(i, curPt);

        //
        // Perform the reflection using the equation
        // R = P-2*N*(dot(N, P)-D)
        // where P = current point, N = planeNormal, D = beta, 
        // and alpha = 2*(dot(N, P)-D)
        //
        double alpha = 0;
        for (int i = 0; i < 3; ++i)
        {
           alpha += planeNormal[i]*curPt[i]; 
        }

        alpha -= beta;
        alpha *= 2.0; 

        for (int i = 0; i < 3; ++i)
        {
            reflection[i] = curPt[i] - (planeNormal[i] * alpha);
        }

        outPts->SetPoint(i, reflection);
    }

    //
    // Add our new points to the output data, and clean
    // up memory. 
    //
    out_ds->SetPoints(outPts);
    outPts->Delete();

    if (gridToPointFilter != NULL)
        gridToPointFilter->Delete();

    return out_ds;
}


// ****************************************************************************
//  Method: avtReflectFilter::AxisReflect
//
//  Purpose:
//      Determines what the mesh type is and calls the appropriate reflect
//      routine.
//
//  Arguments:
//      ds      The dataset to reflect.
//      dim     The dimension to reflect in.
//
//  Returns:    The dataset appropriately reflected.
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2002
//
//  Modifications:
//
//    Hank Childs, Thu Jun 24 12:53:47 PDT 2004
//    Reflect vector data as well.
//
//    Hank Childs, Fri Sep  3 12:10:47 PDT 2010
//    Pass Boolean along to ReflectVectorData function.
//
//    Alister Maguire, Wed Apr 11 09:29:49 PDT 2018
//    Changed name to AxisReflect to distinguish from arbitrary
//    plane reflections. 
//
// ****************************************************************************

vtkDataSet *
avtReflectFilter::AxisReflect(vtkDataSet *ds, int dim)
{
    int dstype = ds->GetDataObjectType();
    vtkDataSet *rv = NULL;
    switch (dstype)
    {
      case VTK_RECTILINEAR_GRID:
        rv = AxisReflectRectilinear((vtkRectilinearGrid *) ds, dim);
        break;

      case VTK_STRUCTURED_GRID:
      case VTK_UNSTRUCTURED_GRID:
      case VTK_POLY_DATA:
        rv = AxisReflectPointSet((vtkPointSet *) ds, dim);
        break;

      default:
        debug1 << "Unknown data type for reflect filter = " << dstype << endl;
        break;
    }

    ReflectVectorData(rv, dim, zeroOutVelocitiesOnBoundary);

    return rv;
}


// ****************************************************************************
//  Method: avtReflectFilter::AxisReflectRectilinear
//
//  Purpose:
//      Reflects a rectilinear mesh.
//
//  Arguments:
//      ds      The dataset to reflect.
//      dim     The dimension to reflect in.
//
//  Returns:    The dataset reflected across val.
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Dec 13 16:41:12 PST 2002   
//    Use NewInstance instead of MakeObject in order to match vtk's new api. 
//    
//    Hank Childs, Wed Aug 11 14:43:44 PDT 2004
//    Create ghost nodes for nodes on a reflection plane.
//
//    Hank Childs, Fri Aug 27 15:25:22 PDT 2004
//    Rename ghost data arrays.
//
//    Hank Childs, Tue Jul  5 09:44:27 PDT 2005
//    Don't produce meshes that are "inside-out"/"degenerate". ['6321]
//
//    Alister Maguire, Wed Apr 11 09:29:49 PDT 2018
//    Changed name to AxisReflectRectilinear to distinguish from arbitrary
//    plane reflections, and set ghostsCreated to true.
//
// ****************************************************************************

vtkDataSet *
avtReflectFilter::AxisReflectRectilinear(vtkRectilinearGrid *ds, int dim)
{
    vtkRectilinearGrid *out = (vtkRectilinearGrid *) ds->NewInstance();
    int nPts = ds->GetNumberOfPoints();
    int dims[3];
    ds->GetDimensions(dims);
    out->GetFieldData()->ShallowCopy(ds->GetFieldData());
    out->SetDimensions(dims);
    out->GetPointData()->CopyAllocate(ds->GetPointData());
    out->GetCellData()->CopyAllocate(ds->GetCellData());
 
    //
    // Reflect across X if appropriate.
    //
    bool flipI = false;
    if (dim & 1)
    {
        vtkDataArray *tmp = ReflectDataArray(ds->GetXCoordinates(), xReflect);
        out->SetXCoordinates(tmp);
        tmp->Delete();
        flipI = true;
    }
    else
        out->SetXCoordinates(ds->GetXCoordinates());

    //
    // Reflect across Y if appropriate.
    //
    bool flipJ = false;
    if (dim & 2)
    {
        vtkDataArray *tmp = ReflectDataArray(ds->GetYCoordinates(), yReflect);
        out->SetYCoordinates(tmp);
        tmp->Delete();
        flipJ = true;
    }
    else
        out->SetYCoordinates(ds->GetYCoordinates());

    //
    // Reflect across Z if appropriate.
    //
    bool flipK = false;
    if (dim & 4)
    {
        vtkDataArray *tmp = ReflectDataArray(ds->GetZCoordinates(), zReflect);
        out->SetZCoordinates(tmp);
        tmp->Delete();
        flipK = true;
    }
    else
        out->SetZCoordinates(ds->GetZCoordinates());


    //
    // Copy over the point data.
    //
    int i, j, k;
    vtkPointData *inPD  = ds->GetPointData();
    vtkPointData *outPD = out->GetPointData();
    for (k = 0 ; k < dims[2] ; k++)
        for (j = 0 ; j < dims[1] ; j++)
            for (i = 0 ; i < dims[0] ; i++)
            {
                int idx = k*dims[1]*dims[0] + j*dims[0] + i;
                int oldI = (flipI ? dims[0]-1-i : i);
                int oldJ = (flipJ ? dims[1]-1-j : j);
                int oldK = (flipK ? dims[2]-1-k : k);
                int oldIdx = oldK*dims[1]*dims[0] + oldJ*dims[0] + oldI;
                outPD->CopyData(inPD, oldIdx, idx);
            }
    vtkCellData *inCD  = ds->GetCellData();
    vtkCellData *outCD = out->GetCellData();

    int iEnd = dims[0]-1;
    if (iEnd <= 0)
        iEnd = 1;
    int jEnd = dims[1]-1;
    if (jEnd <= 0)
        jEnd = 1;
    int kEnd = dims[2]-1;
    if (kEnd <= 0)
        kEnd = 1;
    for (k = 0 ; k < kEnd ; k++)
        for (j = 0 ; j < jEnd ; j++)
            for (i = 0 ; i < iEnd ; i++)
            {
                int idx = k*(dims[1]-1)*(dims[0]-1) + j*(dims[0]-1) + i;
                int oldI = (flipI ? iEnd-1-i : i);
                int oldJ = (flipJ ? jEnd-1-j : j);
                int oldK = (flipK ? kEnd-1-k : k);
                int oldIdx =oldK*(dims[1]-1)*(dims[0]-1)+oldJ*(dims[0]-1)+oldI;
                outCD->CopyData(inCD, oldIdx, idx);
            }
    
    // Figure out which octents are present, meaning that nodes that are on
    // the reflection plane are now interior.
    bool doX, doY, doZ;
    HasNeighbor(dim, doX, doY, doZ);

    vtkUnsignedCharArray *gn = vtkUnsignedCharArray::New();
    gn->SetName("avtGhostNodes");
    gn->SetNumberOfTuples(nPts);
    unsigned char *gnp = gn->GetPointer(0);

    //
    // Start off by turning everything to 0.
    //
    for (i = 0 ; i < nPts ; i++)
        gnp[i] = 0;

    //
    // If the min-face is on the axis of reflection, then turn all the nodes
    // on that face to ghost nodes.
    //
    if (doX && Equal(out->GetXCoordinates()->GetTuple1(0), xReflect))
    {
        for (int k = 0 ; k < dims[2] ; k++)
            for (int j = 0 ; j < dims[1] ; j++)
            {
                int idx = 0 + j*dims[0] + k*dims[0]*dims[1];
                avtGhostData::AddGhostNodeType(gnp[idx], DUPLICATED_NODE);
            }
    }
    if (doX && Equal(out->GetXCoordinates()->GetTuple1(dims[0]-1), xReflect))
    {
        for (int k = 0 ; k < dims[2] ; k++)
            for (int j = 0 ; j < dims[1] ; j++)
            {
                int idx = dims[0]-1 + j*dims[0] + k*dims[0]*dims[1];
                avtGhostData::AddGhostNodeType(gnp[idx], DUPLICATED_NODE);
            }
    }
    if (doY && Equal(out->GetYCoordinates()->GetTuple1(0), yReflect))
    {
        for (int k = 0 ; k < dims[2] ; k++)
            for (int i = 0 ; i < dims[0] ; i++)
            {
                int idx = i + 0*dims[0] + k*dims[0]*dims[1];
                avtGhostData::AddGhostNodeType(gnp[idx], DUPLICATED_NODE);
            }
    }
    if (doY && Equal(out->GetYCoordinates()->GetTuple1(dims[1]-1), yReflect))
    {
        for (int k = 0 ; k < dims[2] ; k++)
            for (int i = 0 ; i < dims[0] ; i++)
            {
                int idx = i + (dims[1]-1)*dims[0] + k*dims[0]*dims[1];
                avtGhostData::AddGhostNodeType(gnp[idx], DUPLICATED_NODE);
            }
    }
    if (doZ && Equal(out->GetZCoordinates()->GetTuple1(0), zReflect))
    {
        for (int j = 0 ; j < dims[1] ; j++)
            for (int i = 0 ; i < dims[0] ; i++)
            {
                int idx = i + j*dims[0] + 0*dims[0]*dims[1];
                avtGhostData::AddGhostNodeType(gnp[idx], DUPLICATED_NODE);
            }
    }
    if (doZ && Equal(out->GetZCoordinates()->GetTuple1(dims[2]-1), zReflect))
    {
        for (int j = 0 ; j < dims[1] ; j++)
            for (int i = 0 ; i < dims[0] ; i++)
            {
                int idx = i + j*dims[0] + (dims[2]-1)*dims[0]*dims[1];
                avtGhostData::AddGhostNodeType(gnp[idx], DUPLICATED_NODE);
            }
    }

    out->GetPointData()->AddArray(gn);
    gn->Delete();
    ghostsCreated = true;

    return out;
}


// ****************************************************************************
//  Method: avtReflectFilter::ReflectDataArray
//
//  Purpose:
//      Reflects a single data array across a line.  This is used as a
//      convenience function for reflecting rectilinear grids.
//
//  Arguments:
//      coords  The coordinates to reflect.
//      val     The point to reflect about.
//
//  Returns:    The coordinates reflected; this must be deleted by the calling
//              object.
//
//  Programmer: Hank Childs
//  Creation:   March 8, 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Dec 13 16:41:12 PST 2002   
//    Use NewInstance instead of MakeObject in order to match vtk's new api. 
//    
//    Hank Childs, Fri Jul  1 10:50:58 PDT 2005
//    Reverse the coordinate array as we go.  Otherwise the grid will be
//    inverted.
//
// ****************************************************************************

vtkDataArray *
avtReflectFilter::ReflectDataArray(vtkDataArray *coords, double val)
{
    //
    // Reflect the coordinates across the line of reflection.
    //
    vtkDataArray *newcoords = coords->NewInstance();
    int nc = coords->GetNumberOfTuples();
    newcoords->SetNumberOfTuples(nc);

    if (coords->GetDataType() == VTK_FLOAT)
    {
        float *c = (float*)coords->GetVoidPointer(0);
        float *n = (float*)newcoords->GetVoidPointer(0);
        for (int i = 0 ; i < nc ; i++)
        {
            n[nc-i-1] = 2*val - c[i];
        }
    }
    else if (coords->GetDataType() == VTK_DOUBLE)
    {
        double *c = (double*)coords->GetVoidPointer(0);
        double *n = (double*)newcoords->GetVoidPointer(0);
        for (int i = 0 ; i < nc ; i++)
        {
            n[nc-i-1] = 2*val - c[i];
        }
    }


    return newcoords;
}


// ****************************************************************************
//  Method: avtReflectFilter::AxisReflectPointSet
//
//  Purpose:
//      Reflects a dataset that is derived from point set.
//
//  Arguments:
//      ds      The dataset to reflect.
//      dim     The dimension to reflect in.
//
//  Returns:    The dataset reflected across val.
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Dec 13 16:41:12 PST 2002   
//    Use NewInstance instead of MakeObject in order to match vtk's new api. 
//    
//    Hank Childs, Wed Aug 11 14:43:44 PDT 2004
//    Create ghost nodes for nodes on a reflection plane.
//
//    Hank Childs, Fri Aug 27 15:25:22 PDT 2004
//    Rename ghost data arrays.
//
//    Brad Whitlock, Wed Apr 11 23:29:56 PDT 2012
//    Preserve coordinate type.
//
//    Alister Maguire, Tue Nov 15 11:51:26 PST 2016
//    Added DeepCopy of ds if threads are greater
//    than 1. 
//
//    Alister Maguire, Wed Apr 11 09:29:49 PDT 2018
//    Changed name to AxisReflectPointSet to distinguish from arbitrary
//    plane reflections, and set ghostsCreated to true. 
//   
// ****************************************************************************

vtkDataSet *
avtReflectFilter::AxisReflectPointSet(vtkPointSet *ds, int dim)
{
    vtkPointSet *out = (vtkPointSet *) ds->NewInstance();

    if (VisItInit::GetNumberOfThreads() <= 1)
        out->ShallowCopy(ds);
    else
        out->DeepCopy(ds);

    vtkPoints *inPts  = ds->GetPoints();
    vtkPoints *outPts = vtkPoints::New(inPts->GetDataType());
    int nPts = inPts->GetNumberOfPoints();
    outPts->SetNumberOfPoints(nPts);

    // Figure out which octents are present, meaning that nodes that are on
    // the reflection plane are now interior.
    bool doX, doY, doZ;
    HasNeighbor(dim, doX, doY, doZ);

    vtkUnsignedCharArray *gn = vtkUnsignedCharArray::New();
    gn->SetName("avtGhostNodes");
    gn->SetNumberOfTuples(nPts);
    unsigned char *gnp = gn->GetPointer(0);

    for (int i = 0 ; i < nPts ; i++)
    {
        double pt[3];
        inPts->GetPoint(i, pt);
        if (dim & 1)
        {
            pt[0] = 2*xReflect - pt[0];
        }
        if (dim & 2)
        {
            pt[1] = 2*yReflect - pt[1];
        }
        if (dim & 4)
        {
            pt[2] = 2*zReflect - pt[2];
        }
        outPts->SetPoint(i, pt);

        gnp[i] = 0;
        if (doX && Equal(pt[0],xReflect))
            avtGhostData::AddGhostNodeType(gnp[i], DUPLICATED_NODE);
        if (doY && Equal(pt[1],yReflect))
            avtGhostData::AddGhostNodeType(gnp[i], DUPLICATED_NODE);
        if (doZ && Equal(pt[2],zReflect))
            avtGhostData::AddGhostNodeType(gnp[i], DUPLICATED_NODE);
    }

    out->SetPoints(outPts);
    outPts->Delete();
    out->GetPointData()->AddArray(gn);
    gn->Delete();
    ghostsCreated = true;

    return out;
}


// ****************************************************************************
//  Function: ReflectVectorData
//
//  Purpose:
//      Reflects the vector data.
//
//  Programmer: Hank Childs
//  Creation:   June 24, 2004
//
//  Modifications:
//
//    Hank Childs, Fri Sep  3 12:10:47 PDT 2010
//    Add support for zero-ing out velocity field along the boundary.
//    This will make it be a continuous field (and one that particles
//    won't advect through!).
//
//    Brad Whitlock, Wed Apr 11 23:28:52 PDT 2012
//    Switch to Get/SetComponent.
//
// ****************************************************************************

static void
ReflectVectorData(vtkDataSet *ds, int dim, bool zeroOutVelocitiesOnBoundary)
{
    if (dim == 0)
    {
        // No reflection.
        return;
    }

    // First identify which variables are vectors.  Point vectors only.
    std::vector<vtkDataArray *> pt_vectors;
    vtkPointData *pd = ds->GetPointData();
    for (int i = 0 ; i < pd->GetNumberOfArrays() ; i++)
    {
        vtkDataArray *da = pd->GetArray(i);
        if (da->GetNumberOfComponents() == 3)
        {
            pt_vectors.push_back(da);
        }
    }
    
    // Now reverse the orientation of the vectors based on the axis of
    // reflection.  Point vectors only.
    std::vector<vtkDataArray *> new_pt_vectors;
    for (size_t i = 0 ; i < pt_vectors.size() ; i++)
    {
        vtkDataArray *da = pt_vectors[i]->NewInstance();
        da->DeepCopy(pt_vectors[i]);
        da->SetName(pt_vectors[i]->GetName());
        const size_t ntups = da->GetNumberOfTuples();
        for (size_t j = 0 ; j < ntups ; j++)
        {
            if (dim & 1)
                da->SetComponent(j, 0, -da->GetComponent(j, 0));
            if (dim & 2)
                da->SetComponent(j, 1, -da->GetComponent(j, 1));
            if (dim & 4)
                da->SetComponent(j, 2, -da->GetComponent(j, 2));
        }
        new_pt_vectors.push_back(da);
    }

    //
    // If we want a continuous velocity field (probably for streamlines), then zero out
    // the velocities on the boundary.
    //
    if (zeroOutVelocitiesOnBoundary && (new_pt_vectors.size() > 0))
    {
        vtkUnsignedCharArray *gn = (vtkUnsignedCharArray *) pd->GetArray("avtGhostNodes");
        double zeroVel[3] = { 0,0,0 };
        size_t npts = ds->GetNumberOfPoints();
        for (size_t i = 0 ; i < npts ; i++)
        {
            if (gn->GetValue(i) != 0)
            {
                for (size_t j = 0 ; j < new_pt_vectors.size() ; j++)
                    new_pt_vectors[j]->SetTuple(i, zeroVel);
            }
        }
    }

    // Now replace all of the original vectors with our new ones.
    for (size_t i = 0 ; i < new_pt_vectors.size() ; i++)
    {
        bool isActiveVector = false;
        vtkDataArray *da = pd->GetVectors();
        if (da != NULL)
           if (strcmp(da->GetName(), new_pt_vectors[i]->GetName()) == 0)
               isActiveVector = true;
        pd->RemoveArray(new_pt_vectors[i]->GetName());
        if (isActiveVector)
            pd->SetVectors(new_pt_vectors[i]);
        else
            pd->AddArray(new_pt_vectors[i]);
        new_pt_vectors[i]->Delete();
    }

    // First identify which variables are vectors.  Now the cell vectors
    std::vector<vtkDataArray *> cell_vectors;
    vtkCellData *cd = ds->GetCellData();
    for (int i = 0 ; i < cd->GetNumberOfArrays() ; i++)
    {
        vtkDataArray *da = cd->GetArray(i);
        if (da->GetNumberOfComponents() == 3)
            cell_vectors.push_back(da);
    }
    
    // Now reverse the orientation of the vectors based on the axis of
    // reflection.  Cell vectors only.
    std::vector<vtkDataArray *> new_cell_vectors;
    for (size_t i = 0 ; i < cell_vectors.size() ; i++)
    {
        vtkDataArray *da = cell_vectors[i]->NewInstance();
        da->DeepCopy(cell_vectors[i]);
        da->SetName(cell_vectors[i]->GetName());
        const size_t ntups = da->GetNumberOfTuples();
        for (size_t j = 0 ; j < ntups ; j++)
        {
            if (dim & 1)
                da->SetComponent(j, 0, -da->GetComponent(j, 0));
            if (dim & 2)
                da->SetComponent(j, 1, -da->GetComponent(j, 1));
            if (dim & 4)
                da->SetComponent(j, 2, -da->GetComponent(j, 2));
        }
        new_cell_vectors.push_back(da);
    }

    //
    // If we want a continuous velocity field (probably for streamlines), then zero out
    // the velocities on the boundary.
    //
    if (zeroOutVelocitiesOnBoundary && (cell_vectors.size() > 0))
    {
        debug1 << "Zero-ing out vectors along the boundary for cell-centered "
               << "vectors.  We have a slow implementation.  If you are looking "
               << "for performance issues, this is the place to look." << endl;
        vtkUnsignedCharArray *gn = (vtkUnsignedCharArray *) pd->GetArray("avtGhostNodes");
        double zeroVel[3] = { 0,0,0 };
        size_t ncells = ds->GetNumberOfCells();
        for (size_t i = 0 ; i < ncells ; i++)
        {
            vtkCell *cell = ds->GetCell(i);
            vtkIdList *ids = cell->GetPointIds();
            size_t nids = ids->GetNumberOfIds();
            bool touchesBoundary = false;
            for (size_t j = 0 ; j < nids ; j++)
            {
                if (gn->GetValue(ids->GetId(j)) != 0)
                   touchesBoundary = true;
            }
            // inefficient
            if (touchesBoundary)
            {
                for (size_t j = 0 ; j < new_cell_vectors.size() ; j++)
                    new_cell_vectors[j]->SetTuple(i, zeroVel);
            }
        }
    }

    // Now replace all of the original vectors with our new ones.
    for (size_t i = 0 ; i < new_cell_vectors.size() ; i++)
    {
        bool isActiveVector = false;
        vtkDataArray *da = cd->GetVectors();
        if (da != NULL)
           if (strcmp(da->GetName(), new_cell_vectors[i]->GetName()) == 0)
               isActiveVector = true;
        cd->RemoveArray(new_cell_vectors[i]->GetName());
        if (isActiveVector)
            cd->SetVectors(new_cell_vectors[i]);
        else
            cd->AddArray(new_cell_vectors[i]);
        new_cell_vectors[i]->Delete();
    }
}


// ****************************************************************************
//  Method: avtReflectFilter::HasNeighbor
//
//  Purpose:
//      Determines if this reflected domain has a neighbor across each of the
//      three reflection planes.  If it does, then it will be internal to the
//      problem.  If not, it will still be external.
//
//  Programmer: Hank Childs
//  Creation:   August 11, 2004
//
// ****************************************************************************

void
avtReflectFilter::HasNeighbor(int octent, bool &doX, bool &doY, bool &doZ)
{
    const int *shouldreflect = atts.GetReflections();

    int xNeighbor = ((octent & 1) ? octent-1 : octent+1);
    doX = (shouldreflect[xNeighbor] != 0 ? true : false);

    int yNeighbor = ((octent & 2) ? octent-2 : octent+2);
    doY = (shouldreflect[yNeighbor] != 0 ? true : false);

    int zNeighbor = ((octent & 4) ? octent-4 : octent+4);
    doZ = (shouldreflect[zNeighbor] != 0 ? true : false);
}


